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ABSTRACT 

One of the most intriguing scenarios proposed to explain how active galactic 
nuclei are triggered involves the existence of a supermassive binary black hole 
system in their cores. Here we present an observational evidence for the first 
spectroscopically resolved sub-parsec orbit of a such system in the core of Seyfert 
galaxy NGC 4151. Using a method similar to those typically applied for spectro- 
scopic binary stars we obtained radial velocity curves of the supermassive binary 
system, from which we calculated orbital elements and made estimates about the 
masses of components. Our analysis shows that periodic variations in the light 
and radial velocity curves can be accounted for an eccentric, sub-parsec Keple- 
rian orbit of a 15.9-year period. The flux maximum in the lightcurve correspond 
to the approaching phase of a secondary component towards the observer. Ac- 
cording to the obtained results we speculate that the periodic variations in the 
observed Ha line shape and flux are due to shock waves generated by the su- 
personic motion of the components through the surrounding medium. Given the 
large observational effort needed to reveal this spectroscopically resolved binary 
orbital motion we suggest that many such systems may exist in similar objects 
even if they are hard to find. Detecting more of them will provide us with insight 
into black hole mass growth process. 

Subject headings: galaxies: active — galaxies: interactions — (galaxies:) quasars: 
individual (NGC 4151) — shock waves — galaxies: Seyfert — black hole physics 
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Introduction 



Dif 
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2003 




Gaskell 


2009 
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for e.g. 



cresting possib i 



Komossa 



Eracleous et al. 



2006 



2012 



and the references therein). One of the 
ities involves the existence of a binary system in the core (see 



Bog danovic et al. 



Popovic 



2012 



2008 



Gaskell 



2009 



Tsalmantza et al. 



2011 



and the references therein). If black hole (BH) mass 
grows via major mergers we might expect to see the signature of a binary black hole in 
some or many active galaxies, if the merger process involves slow coalescence of the two 
components. 

NGC 4151, one of the best studied Seyfert galaxies, shows a complex sub-parsec 
structure. Studies at radio, visible, and X-ray wavelengths indicate violent processes 
including ejection of gas. Detection of outflow scales ranging from several hundred 



Kraemer et al. 


2001; 


Mundell et al. 


2003; 


Kraemer et al. 


2008; 


Mundell et al. 


1999 



and 



the references therein). Vio lent variability in this act ive galactic nucleus (AGN) has been 



monitored over many y ears (S hap ova 



ova et al. 



20081 ). Earlier work showed variation in the 



wings of the CIV line ( jUlrich et al. Ill985l ) that were associated to two localized emitting 
regions with line of sight velocities of -6100 and +8500 km s _1 with respect to the systemic 
velocity. These localized regions were thought to be the possible signatures of orbiting 



matter or of a two-sided jet. A radio jet or highly inclined disk 



perpendicular to the arcsec-scale radio jet was imaged at 6 cm (jUlvestad et al 



length ~ 1 pc) 



align ed 



1998|). 



Flux variations in the optical domain have been monitored over a m uch longer time (over 



100 years) (see for e.g. 



Guo et al. 



2006 



Oknyanskij fc Lyuty 



2007, and the references 



therein). Optical and UV emission is unresolved so clues about central structure come 
mainly from analysis of broad emission line profile. Reverberation mapping studies (see 
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for e.g. lOnken et al. 1 120071 ; IShapovalova et al. 1120081 and the references therein) suggest a 
very compact broad line region (BLR) with broad Balmer lines responding to continuum 
changes within a few days. Many authors claim to find long term p eriodicities in the 



flux variations of NGC 4151 (jLongo et al. 



1978 



Oknyanskij fc Lyuty 



2007 



Guo et al. 



996) of e.g. 15.8 years 



2006 



Pacholczyk et al. 



Okn yanskij et al 



1983 



Chuvaev et al. 



20081 ). Compared to Seyfert Is, NGC 4151 has a harder X- 



black-hole binaries in the hard state (ILubinski et al. 



2010|). 



ray spectrum and also similar to 



In this paper we analyze the Ha line shape and flux variability in NGC 4151 during 
period of more than 20 years in order to investigate the BLR structure in NGC 4151. The 
possibility of periodicity in the line profile variations leads us to interpret it in terms of 
binary system orbital motion at sub-parsec scale in the center of NGC 4151. 

The paper is organized as follows: in Section 2 we briefly describe our data set, 
method, and analysis used to detect periodicity in both the line flux and the shape of the 
profiles. In Section 3, we discus possible scenarios for the obtained periodic variations, and 
give evidence of a sub-parsec supermassive binary BH system in the core of NGC 4151. In 
Section 4, we discuss alternative models, and implications of our results. Finally, in Section 
5, we summarize results and give conclusions. 



2. Analysis 



We analyzed 115 spectra of the Ha line profile in NGC 4151 covering more than 20 
years. Most of the spectra were obtained during a period of 11 years from 1996 to 2006 



( IShapovalova et al. 



2008 



20101 ). This dataset was supplemente d with: (i) AGN w atch 



spectra observed for three months beginning 1993 November 14 (IKaspi et al. 



199J) with 
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the lm Wise Telescope Q, (ii) one spectrum from Asiago 1.82 m Ekar telescope of the 
Padova Astronomical Observatory (1989 March 24), which was equipped with a Boiler & 
Chivens spectrograph and a 600 gr/mm grating yielded a resolution of 4- 5 A FWHM, wi th 



1995 ^ 



PA=90, at the Cassegrain focus and (iii) a spectrum from 1986 March 29 (IHo et al. 
All spectra were reduced to a common resolution of 15 A (about a = 340 km s -1 ). The 
continuum was removed by subtracting a linear fit between 6265 and 6830 A. All spectra 
were scaled to the flux of [01] 6300 A which is assumed constant at all epochs. 



2.1. Flux Variations 



We derived Ha light curves for the total flux (see Fig. 1) as well as for the partial flux 



in 1000 km s 1 ve locity bins ( 



spectral analysis (ILomb 



1976 



ig. 2). We an alyzed the light curves using Lomb-Scargle 



Scargle 



19821 ) procedures. 



The Lomb-Scargle normalized periodogram (ILomb 



Lomb 


1976; 


Scarele 


1982; 


Press et al. 



19961 ) is a powerful tool for finding and testing the significance of weak periodic signals 
embedded in otherwise random and unevenly sampled observations. The significance level 
(i.e. the probability that the data contain a periodicity) can be tested fairly rigorously by 
the statistic 1 — p, where p is the so-called "false alarm probability" . This estimates the 
probability that the peak occurred in the presence of pure independently and normally 
(Gaussian) distributed noise. A small value for the false-alarm probability indicates a 
significant periodic signal. Resultant periodicities were subsequently tested by fitting with 



x The AGN watch data could be obtained in the digital format from the following link: 



http:/ /www. astronomy.ohio-state.edu/ ~agnwatch/data.html 



2 The spectrum from Ho et al. 1995 is public and could be obtained in digital format from 
the NED: http://ned.ipac.caltech.edu 
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sine functions (see Fig. 1, bottom panel). 

Since, we detected the periodicity in the Ha line lightcurve (see Fig. 1), we tried to 
identify the parts of the line profiles that are most affected by those periodical variations. 
For that purpose we analyzed the light curves for each 1000 km s _1 radial velocity bin of 
Ha spectra, with significance level above 90 % (see Fig. 2). Since, we detected periodicity 
in certain light curves (mainly in the core of the line for the velocity bins between -3000 
and 4000 km s _1 ), we tried to uncover possible mechanisms that could produce such flux 
variations. For that purpose we tested many different line decomposition models, taking 
into account the parts of the line where we detected significant periodicity. 



2.2. Radial Velocity Curves 

We used Gaussian analysis to measure profile shape variations. This was accomplished 
by fitting each Ha profile with the sum of Gaussians assuming some constrains. The 
intensity ratio of narrow components (Ha line, [OI]6300,6364 A, [NII]6548,6584 A, 
[SII]6717,6731 A and Hel 6678 A ) were modeled using spectra obtained during the 
minimum activity stage (2005 May 12). They were fitted with sum of Gaussian profiles 
assuming that the intensity ratios of the narrow lines did not change during the monitoring 
period. The widths of the narrow lines were assumed equal during each fit. 



The broad Ha profile was fit together with the narrow lines. We tested a number of 



mode ls including those with more than one Gaussian component (see for e.g. 



Shen&Loeb 



20101 ) and analyzed the light curves of the component fluxes as well as the radial velocity 



curves of each component during the monitoring period. The choic e of component width s 



was not arbitra r y. W e started from the decomposition proposed i nlSulentic et al. 



Marziani et al 



(120101 ). We noticed that a two component model (IBon et al. 



2009 



(2000); 



200J) 



-7- 



consisting of a very broad (covering the far wings) and broad (covering the line core) 
components could not properly explain an intermediate width feature visible on the Ha 



wing (a = 600 


mi s" 


(Sulentic et al. 


2000 



tm s 1 ). R and VB 



Marziani et al 



observed in this class (jZamfir et al. 



R co mponents are typical of Population B AGN 



2010) while displaced subpeaks have been occasionally 
201ol ). The intermediate width subpeak was observed 
on both the blue and red wings of Ha at different epochs but never at the same time leading 
us to interpret it as a single emitting component (a = 600 km s^ 1 ). 

We then tried many fits with different widths for each of the three Ha components 
and adopted a model that consists of: (i) a very broad component (VBC) corresponding to 
the far wings (a = 3400 km s _1 ), (ii) a central broad component (CBC) corresponding to 
the "classical" broad core of the line (a = 1700 km s _1 ) and (iii) an intermediate width 
component (bump) appearing on either the blue or red wing of the Ha profile a = 600 
km s _1 ). The choice of Gaussian component widths were made in such a way that each 
component width differs for at least two times, providing a clear identification of each 
component. We fit the Ha profile using three components with fixed widths. The CBC 
shift was fixed to the shift of the narrow lines. 

To secure the identification of the Gaussian components we adopted the model with 
three Gaussians with several important limitations according to the shapes of profiles. We 
measured the width of the bump on profiles where it could be clearly isolated from other 
components. We found that the width of the bump's Gaussian is around a=600 km s _1 . 
We also tested the widths of other two components and found that the Gaussian suitable 
to cover the central part of the line profiles (CBC) needed to be at least twice the width of 
the bump's Gaussian. Also, we fixed the shift of the CBC since we noticed that the shift 
of this component was showing very small variations (smaller than the estimated errors 
for this parameter), since the very broad wings (which were significantly building up in 
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many profiles) could be covered only with very broad Gaussian (at least twice the width of 
the CBC). We noticed that when the bump component is present on the red side of the 
profile, it was appearing either as a separated peak, or as feature that was making the red 
side of the CBC much steeper than the regular shape of the Gaussian component that fit 
the core of the line (CBC). When the bump component was appearing on the blue side, it 
was just changing the steepness of the blue side of the profile, making it much steeper than 
the shape of the Gaussian of the CBC. In this way including the bump component was 
justified in our model. Similarly, it was necessary to introduce the VBC Gaussian to fit the 
far wings in some profiles, that also could not be reproduced by other two components. 

In the case where we vary the intensities of the CBC and the VBC, the position of the 
bump would shift, but for the solutions where the fit looks reasonable, the shifts of the 
bump are smaller than the estimated errors. For example, if intensities are changed for less 
than 5%, the shift of the bump is changed, but stays inside the error interval, while reduced 
X 2 increases for more than an order of magnitude. If the variations of intensities are even 
larger the fit becomes unreasonable, and do not follow the shape of the profile. 

Examples of fits for two epochs are shown in Fig. 3. We noticed that the CBC and 
VBC give major contributions to the total flux of the line, and hence their lightcurves follow 
the total Ha line lightcurve. On these lightcurves (see Fig. 4. bottom) one can identify one 
dominant peak in the interval of the Modified Julian Date (MJD) between 50000 and 51000 
MJD and another smaller peak between 52500 and 53000 MJD. Since we found similar 
behavior in the lightcurves of all components, as well as in the continuum and the Ha line, 
we concluded that the flux variations are real in all components, and are not the caused 
by of the Gaussian fitting procedure. More spectra from each observational campaign with 
corresponding fits are presented in Figures 5-10 



Bogdanovic et al. 



(120081 ) made simulations of Ha emission line profiles in the case of 
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a supermassive binary black hole system with masses of 10 7 and 10 8 M Q and the following 
orbital elements: 15.7 year period and semimajor axis of a = 0.01 pc. It is not that obvious 
but one can see in their Figure 18 that the emission line profiles show bump peak that 
appears on the blue side of the profile for the first part of the orbital period, and in the 
second part of the orbital period the bump peak starts to appear on the red side of the 
profile. Also, one can see that the central peak is present in the most of profiles as well 
as the wider structure (VBC) even though we approximated it with a Gaussian instead of 
more complex profile as seen in this figure (probably a disk like emission profile). In our 
case there was no significant difference in the radial velocity curves and light curves in the 
case of fitting using a VBC Gaussian instead of a VBC disk profile. 

We analyzed radial velocity variations in the VBC and bump compone nts. We also fi t 



radial velocity curves assuming Keplerian orbits using the "Velocity" code (jRainer II1988I ). 
The uncertainty in the position of the bump was estimated to be ± 200 km s _1 (when 
detected on the blue side) and ± 70 km s" 1 (on the red) at 2a confidence level (P=0.99). 
Estimated uncertainties for radial velocity of the VBC were around ± 300 km s^ 1 . These 
errors were set by x 2 behavior when varying only the displacement from the best fit radial 
velocity and are regarded as minimum uncertainties. 

To test the influence of the fitting procedure on the lightcurves and radial velocity 
curves we varied the widths of each component. For the bump we found that in some 
cases the fit becomes unreasonable if we change its width by more than 100 km s _1 . For 
the change in the widths of the VBC and CBC, the radial velocity curves of the bump 
component did not change much (a few times less than the estimated errors), even thought 
the VBC width was varied from 3000 to 4500 km s -1 , and CBC from 1500 to 2000 km 
s _1 . We noticed that for the broader VBC, the amplitude of the radial velocity curve of 
that component, became larger, but kept the same shape and behavior, and hence led to 
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the similar orbital parameters. However, in that case the obtained mass ratio between 
secondary (bump) and primary (VBC) will vary between 0.2 (for ovbc = 3000 km s _1 ) and 
0.6 (for cr V Bc = 4500 km s _1 ). Therefore, this method gives more reliable estimations of the 
orbital elements than those for the mass ratio. 

We are aware that, when the bump is close to the line core (at v r ~ 0) its position 
is well constrained because of the model we build with three symmetric Gaussians, since 
the bump is overwhelmed by the BC and VBC emission so as to become not visually 
identifiable. As the bump radial velocity changes sign rather suddenly, this condition occurs 
only for a few epochs. 



3. Results 

Light curve analysis is presented in Figs. 1, and 2. Fig. 1 shows that the total flux 
could be well fit with a sine function with 5700 ± 300 day periodicity. Fig. 2 shows 
Lomb-Scargle analysis fluxes for each 1000 km s _1 radial velocity bin of Ha with significance 
level above 90 %. The resultant values ranged from P « 5000 to 6000 days (see Fig. 2). 



3.1. Periodical Variations in Lightcurves and Possible Mechanisms 



The detection of a significant periodicity is a major result of this study. The analysis 
of the light curves showed a periodicity of nearly 16 years, which is in agreement with 
previous studies (of 16 year periodicity) obtained (using different methods) for a much 



earlier monitoring period, covering few earlier cycles (see 



Oknyanskij et al 



19781 ). 



There are several possible sources of periodicity. Apart from the orbital period 

_ i _ i 

(Tp m 5.5r 16 2 M 7 2 yr), relevant timescales include the timescale associated with geodetic 
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precession Tp w 2 ■ 10 4 r 16 2 (— )M 7 2 yr (i.e., precession o f the BH spin around the total 



angular momentum of the system (IBegelman et al 



1980)), and accretion disk precession 



due to a second massive b 



the accretion disk radius (IKatz 



i — - — - 

ack hole T P A n ~ 23 ■ rf 6 (^) 2m 6 2 R d f 6 cos" 1 6> yr, where is 



19971 ). Measured parameters for NGC 4151 imply that 
the geodetic precession is Tp ~ 10 5 yr. Accretion disk precession is expected in any AGN 
hosting a binary black hole and Tp^D is much lower than the value estimated for geodetic 

_ 3 

precession. However for NGC 4151 Tp^ad ~ 250f? d ^ 6 cos" 1 9q is still significantly longer than 
the derived period. A third possibility is that radiation force induce a "self-warping" of 



the accretion disk 



would be (IPringle 



(obviating the necessity for a second black hole). In this case the period 

i i 

19961) T PiWarp « 13.6 • (M ion . gas _ 2 )M 7 2 i? 2 10 3 R L± 4 yr. The self-warping 



timescale is consistent with the observed P and could, in principle, reinforce a periodic 
behavior also involving continuum fluxes. 



3.2. A Sub-parsec Supermassive Black Hole Binary 

We analyzed radial velocity curves for the VBC and bump components, because we 
noticed large variations of their velocity shifts. Fig. 4 indicates that the radial velocity 
curve can be well described with Keplerian orbits in a binary system. 

The derived orbit is eccentric (e = 0.42) with a period of P = 5776 days and longitude 
of pericenter to « 95°. We note that orbital precession would be expected but we did not 
take it into account since the predicted 5$ < 1° per period in Schwarzschild metric. Using 
Kepler's laws we estimated the semimajor axes: aisinz = 0.002 pc, a2sinz = 0.008 pc 
and masses: mi sin 3 i — 3 • 1O 7 M and rri2 sin 3 i — 8.5 • 1O 6 M (where i is the inclination 
of the orbital plane). These results are in a good agree ment with the correspon ding; 



theoretical results obtained using numerical simulations ( IBogdanovic et al. 



20081 ) which 



showed simulations (including the emission line profiles) of an eccentric Keplerian system 
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with 15.7 year period, at 0.01 pc distance with similar total mass as found here. We 
also, estimated semimajor axes and masses for different inclination angles and found that 
an inclination of about 65° yields semimajor axes of 0.0024 pc and 0.0094 pc, while the 
masses of each component are 4.4 • 10 7 M Q and 1.2 • 10 7 M & . Our mass estimate is in 
good agreement wit h the values obtained from other methods e.g. reverberation mapping 



( lOnken et al. 



20071 ). The derived separation of two supermassive black holes (SMBHs) 



of about 0.01 pc is expected fr om theoretical models for the evolution of binary SMBHs 



( IMilosavljevic &: Merritt 



20011 ). Th e derived period is in agreement with previous stud ies 



of optical variab ility in NGC 4151 ( jOknyanskij et al 



1978 



Guo et al. 



Oknvanskii fe Lvutv 



20061 ) . which indicate existence of previous cycles (see 



2007 



Oknyanskij et al 



19781 ). 



If this model is correct then gravitational radiation should lead to coalescence 



( jMisner et al. 



1973 ) in ~ 10 8 yr. The evolution of a BBE 



topic of much recent research (see 



Merritt &: Milosavljevic 



syst 



em through merging is a 



20051 . for an exhaustive review) 



Regardless of the evolutionary stages governed by stellar dynamics, there will be a final 
stage when dissipation of orbital energy is due to emission of gravitational radiation. 
When the binary separation is on sub-parsec scales, coalescence time scales oc M~ 3 imply 
that such binary configurations co uld be relatively long- lived in low- mass AGNs, if other 



mechanisms do not intervene (e.g. 



Hayasaki et al. 



20081 ). Systems like NGC 4151 with a 



central supermassive black hole and a (or perhaps, even more than one) much less massive 
satellite could be common even if they are hard to find. 

Considering the relatively large uncertainty in radial velocity for the VBC, the mass 
ratio could be different than we have derived. In the most extreme case we could be 



observing some cloud of gas (jGillessen 



20121 ) spiraling down toward the center of attraction, 



but then we should not expect a very long lifetime for such cloud due to dissipation 
making this scenario less probable, since several 16 year cycles are already reported 
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( jOknyanskij et al 



1978 



Oknyanskij fe Lyuty 



2007 



Guo et al. 



2006 



including the result 



of this paper). If the secondary object is compact but with a mass much smaller than 
the primary component, then its orbit would be defined with the bump component only, 
implying a distance of about 0.01 pc and the same period (about 15.9 years and eccentricity 
of about 0.4). 



4. Discussion 



If we assume that a second massive black hole is i ndeed present, then the pre tty large 



value of eccentricity suggests an analogy with OJ 287 (IValtonen fc Ciprini 



20121 1 . where 



its orbital motion is assumed to be an explanat ion of its flux varia tion. Orbital motion of 



SMBHs has be en recently reported in 3C 66B (ISudou et al. 



(jGaskell 



2003|) and earlier in 3C390.3 



19961 ) . We are however, far from the extreme conditions proposed for OJ 287, a 



source with a black hole mass estimated to be ~500 times larger then the one of NGC 4151. 
As a consequence, we expect a small pericenter precession. Energy loss due to gravitational 
radiation is very small if compared to OJ 287. Nevertheless, this system is among the 
ones where further evolution is governed by gravitational losses. Since these systems are 
expected to occur when the binary separation roughly corresponds to the size of the BLR, 
extensive monitoring may reveal other cases (we remark that the observational coverage of 
NGC 4151 is almost unique). 

If the component separation in a binary SMBH system is sub-parsec: i.e. on the same 
order or smaller than the expected radius of the BLR, we have justification for a CBC 
fixed to the rest frame as derived from the redshift of the narrow lines. This implies that 
the orbit of the binary SMBH system is most probably located within the BLR region. 



The estimated orbital velocities (v 



v r ■ sin 1 i ■ sin 



where % represents inclination 



and the orbital phase) are above 4000 km s 1 for the bump and around 1000 km s 1 
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for VBC for the quadrature phase (where sine/) = 1 and the radial component is highest 



towards the observer). These are much larger than the typical spe ed of sound (typica 



order of tens to hundreds of km s exp ected for such conditions (jBogdanovic et al. 



Dopita M. A. 



1995 



Maver et al. 



shock waves in the BLR (jShi et al 



ly o 



2008 



2007) . Therefore, one would 



201J : 



Artymowicz fc Lubow 



1996 



them to proc 


uce 


Maver et al. 


2007) 



The epoch with maximum total flux (see Fig. 1, bottom panel), in the MJD interval 
between 50000 and 51000, corresponds to the epoch when the bump (secondary) component 
is approaching the observer with supersonic velocity, producing shock waves in the direction 
of the observer (see Fig. 4), while a smaller flux maximum that follows afterwards (MJD 
between 52000 and 53000) corresponds to the epoch when the approaching shock wave is 
produced by the VBC. The epochs where the fluxes are at minimum (Fig. 4) correspond to 



the phase when the shock waves are propagating perpendicula r to the line-of-sight direc 



tion 



(conjunction phase). This is consistent with previous results ( jShapovalova et al. 



2008|): 



non-radiatively heated emission is needed to explain an excess of line emission with respect 
to the case of pure photoionization during the epoch of the flux maximum (MJD between 
50000 and 51000). Shocks can provide a suitable heating mechanism. 

In other scenarios, involving self-warping and flares, one would expect circular 
concentric orbits with sinusoidally shaped radial velocities. This differs from the radial 
velocity curves in Fig. 4 that corresponds to the more eccentric Keplerian orbits. Also, 



the simulated light and radial velo city curves o 



fragmented spiral arm presented in 



Lewis et al. 



an el 



iptical disk, single spiral arm or 



( 120100 . could not describe well the shape 



of the light and radial velocity curves observed in NGC4151. 

One of the alternative explanations for peaks or shoulders seen in the broad Ha 
pro files could be the emissio n of an accretion disk. To test such possibility we tried to fit 



the 



Chen &: Halpern 



( 19891) model of accretion disk into broad line profiles from different 



- 15 - 



epochs, in such a way that the peak of the bump would be fitted with one of two peaks from 
the disk profile. We found that it could be possible to fit an accretion disk model in most 
of the profiles, but always with needing the additional components. For example, we could 
fit the disk profile into the core of the line in such way that the bump component would 
be covered with one of the peaks from the disk. In that case, the result for the inclination 
parameter is very small (around 15 degrees), which is no t that likely 
much higher inclination of jets, around 60-70 degrees (see 



that likely 


' exp 


anation due to 


Das et al. 


2005 




Kraemer et al. 



200ll ). The variations in observed profiles could be explained assuming that the emission 
l ines originate in the precessing accretion-ring, or precessing elliptical accretion-ring 



( jStorchi-Bergmann et al. 



19971 ) . with the variations of the width of the inner radius of 
emitting ring, or with the shift of the disk profile (so the moving bump would be fitted 
with one of peaks from the disk profile). In all these cases it was necessary to involve 
additional components (e.g. VBC and CBC). If the wider disk profiles were used to fit into 
the far wings, then the corresponding inclination would be around 45 degrees, covering 
the far wings of the profiles, as in the VBC in our model, resulting in a similar behavior 
to the VBC. This will include the need for the CBC and a bump component and with 
radial velocity curves similar to those in our simplified model with three Gaussians. We 
found cases with a higher red peak than blue, that could not be explained with the cir cular 



19951). 



disk model, but could be explained with elliptical disk model (lEracleous et al. 
For this alternative model, its precession could explain its rising and moving peaks. All 
alternative explanations must be consistent with the obtained periodicity (see §3.1). Some 



other alternative exp 



e. g. 



Newman et al. 



anations could be: orbital motion of a flare or a bright spot (sec for 



19971 ). self- warping accretion disk ( jPringle Ill996[ ). precession of an 



eccentric circular and elliptical accretion disk, circular disk with sp iral arm (see for e. g. 



Lewis et al. 



20101 ) , the precession of a bipolar o utflow (see for e 



20071 ). and an inhomogeneous disk (see for e. g. 



Eracleous et al. 



Crenshaw and Kramer 



19951 ). For the most of 
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these sce narios we should expect more sinusoidal-like behavior of radial velocity c urves (see 



for e. g. 



Newman et al. 



1997 



Storchi-Bergmann et al. 



1997 



Lewis et al. 



20101 ). which 



differs from Fig. 4. For these alternative models, more tests should be made, to prove or 
disprove any of them. 

The radial velocity curve is immediately consistent with orbital motion as we are able 
to fit it with an eccentric orbit. The idea of orbital motion seems the most plausible. It 
is a consequence of the favorable value of the argument of the periastron (i.e., a fortunate 
coincidence) that we are able to see a sudden change in the radial velocity. It is however 
what is expected and observed in many binaries with highly eccentric orbits. 



5. Conclusion 

We have analyzed the shape and flux variations in the Ha light curve for the Seyfert 
galaxy NGC 4151 using observations spanning 20 years. We also studied variations in 
different parts of the line profile. Lomb-Scargle spectral analysis of these light curves 
suggest variations that can be satisfactorily fit with sine functions and indicating periodicity 
of about 5700 days for the total flux variation and around 5800 days for most of the 1000 
km s _1 Ha intervals. Analysis of radial velocity curves for Ha emission line components 
(VBC and bump) reveal evidence for orbital motion consistent with the signature of a 
sub-parsec scale SMBH system with orbital period of about 5780 days. 

In this scenario periodic variations in the optical spectra of NGC 4151 could be 
generated by supersonic motions of the components through an extended broad line region. 
Maxima in the flux variations correspond to phase when a shock wave is generated toward 
the observer. This mechanism is also able to explain the detected outbursts in the optical 
fluxes of emission lines and the continuum and is in agreement with reported shock heating 
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during the period around the light curve maximum (see IShapovalova et al. 1120081 ). as 
well as the similarity of its X-ray spectra to black-hole binary spectra in the hard state 



( ILubinski et al. 



20101 ) 



This interpretation is consistent with some previous theore tical results (IBogdanovic et al 



2008 



Mayer et al. 



2007 



2010; 



Milosavljevic fc Merritt 



20011 ) and would support the idea 



of BH mass growth by major merger rather than by slower accretion processes. If the 
phenomenology outlined by our analysis is indeed general for all similar AGN then an 
important question arises: Is binarity a necessary condition for AGN ignition? Support for 
such an assumption could be related to changes in spectral type (Syl-Syl.5-Sy2) manifested 
by NG C 4151. The activity type of the galaxy changes from Syl to Syl.5 when it is most 



active (|Shapovalova et al. 



(iPenston and Perez 



1984 : 



2008) and to Sy2 a t its deep minimum phase, as it was in 1984 



Lyutyj et al. 



19841 ) . Different spectral types may be possible 



only at a certain phase of orbital motion in the binary SMBH. This opens the question 
whether the different activity types of active galaxies correspond to different orbital 
phases of such systems and whether the binarity is a necessary condition for the activity 
switch-on. Further monitoring is needed to assess the persistence (or disappearance) of the 
orbiting components. If our results would be confirmed by future monitoring campaigns or 
analysis on different sets of data, then NGC 4151 might become an important candidate for 
gravitational wave detection in the future. 



This research is part of projects 176003 "Gravitation and the large scale structure of 
the Universe" and 176001 "Astrophysical spectroscopy of extragalactic objects" supported 
by the Ministry of Education and Science of the Republic of Serbia, contract P08-FQM-4205 
from La Junta de Andalucia and also grants N09-02-01136a and 12-02-01237a supported by 
RFBR. We also thank Samir Salim and Giovanni La Mura for constructive discussions. 
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Fig. 1. — Lomb-Scargle normalized periodograms of the light curves for total flux (top) and 
the corresponding sinusoidal fit (bottom). 
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Fig. 2. — Period (filled circles with error bars) obtained for each 1000 km s _1 light curve of 
Ha line profile and their normalized powers (asterisks). Horizontal straight lines represents 
50%, 90% and 99% normalized power significance levels. 
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Fig. 3. — Examples for the best fit (dashed line) of the Ha line profile (solid line) for two 
epochs: MJD 51203 (up) and 52621 (down). Three Gaussian components: VBC (a = 3400 
km s -1 , shaded with gray), CBC (a = 1700 km s -1 , shaded with line pattern) and the 
bump (a = 600 km s _1 , shaded with black) component, are fitted in the broad profile. The 
dash-dotted line represents the narrow line template. The panel at the bottom of the plots 
represents residuals of the fit. 
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Fig. 4. — Keplerian orbits of a supermassive black hole binary. Top: radial velocity curves 
of the bump (open triangles with y error bars) and VBC (open circles with y error bars) 
obtained from Gaussian decomposition of the broad Ha line, as well as fitted radial velocity 



curves for orbits of both components (solid anc 
curves for NGC41 51 at 512.5 nm, compiled fr om 



dashed lines). Midd le: continu um flux light 



Malkov et al. 



19961) (c rosses ) , 



Sergeev 



Shapovalova et al. 



(119971 ) (plusses). 



Kaspi et al. 



(120081 ) (open squares), as well as at 656.3 nm from 
(119941 ) (filled circles). Bottom: light curves of the relative fluxes (normalized to 
the flux of 01 line) with y errorbars of VBC (doted line), CBC (dashed line with open 
squares), bump component (dashed line) and the Ha line (solid line). Illustration of the 
corresponding orbital phases are presented above the radial velocity curve panel, assuming 
counter clockwise direction of each black hole motion and the direction toward the observer 



below the plots. 
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Fig. 5. — The same as in Fig. |3]but for the spectrum from 
46518.5). 
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Fig. 6. — The same as in Fig. [3j but for the spectrum obtained on Asiago 1.8 m Ekar 
telescope (epoch MJD 47609.529). 
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ig. 7. — The same as in Fig. [3] but for the spectrum obtained during AGN watch program 



Kaspi et al. 



1996|) (epoch MJD 49327.577). 
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ig. 8. — T he same as in Fig. [3] but for the spectrum from the set of 



2008 



Shapovalova et al. 



2010l ) (epoch MJD 51166.656) with the bump component positioned in the blue part 



of the spectrum. 
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Fig. 9. — The same as for Fig. Elbut for the epoch MJD 52016.36, with the bump component 
positioned in the central part of the spectrum. 
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Fig. 10. — The same as for Fig. Elbut for the epoch MJD 52621, with the bump component 
positioned in the red part of the spectrum. 



